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SUMMARY 

A numerical technique has been developed for the analysis of specified 
transonic airfoils or for the design of airfoils having a prescribed pressure 
distribution, including the effect of weak viscous interaction. The method 
uses the full potential equation, a stretched Cartesian coordinate system, and 
the Nash-Macdonald turbulent boundary layer method. Comparisons with experi- 
mental data for typical transonic airfoils show excellent agreement. An ex- 
ample shows the application of the method to design a thick aft-cambered air- 
foil, and the effects of viscous interaction on its performance are discussed. 

INTRODUCTION 

A numerical method for the design or analysis of transonic airfoils should 
not only be accurate but also should be as simple as possible in concept and 
approach. It should use coordinate systems, input variables, and boundary 
condition treatments that can be easily understood by the user. Finally, it 
should be able to handle both shocked and shockless flows and be suitable not 
only for complete design but also for airfoil modification. 

One approach to this problem is the inverse method in which the airfoil 
surface pressure is specified and the airfoil shape subsequently determined. 
Admittedly, this approach requires knowledge of what would be a desirable 
pressure distribution, but this characteristic is probably understood by the 
designer as well as any other. Furthermore, the designer can select a pressure 
distribution that will approximate a desired lift and moment, have a reasonable 
supersonic zone, yield desirable boundary layer properties, and satisfy other 
transonic flow criteria. 

The purpose of this paper is to present results obtained with a numerical 
method that is suitable for the analysis, design, or modification of subsonic 
and transonic airfoils. The method is similar to that of references 1-3, but 
it has been modified to include the effects of weak viscous interaction. 

SYMBOLS 


Cp drag coefficient 

C lift coefficient 

.Li 

pitching moment coefficient 
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pressure coefficient 
chord 


C 

P 
c 

M free stream Mach number 

00 

RN Reynolds number 
w weighting factor 
a angle of attack 

6* boundary layer displacement thickness 

PROBLEM FORMULATION 

To solve the inviscid part of the flowfield, the method uses the exact 
equation for the perturbation potential in Cartesian coordinates. In order to 
avoid at supersonic points difficulties associated with nonalignment of the 
coordinates and the flowfield, a rotated finite difference scheme (ref. 4) is 
used in the solution; and in the actual program the infinite physical plane is 
mapped to a rectangular computational box. The resulting transformed finite 
difference equations are solved iteratively by column relaxation sweeping from 
upstream to downstream. 

In the design mode, the shape of the nose region (typically 6-10% chord) 
is specified and a pressure distribution is prescribed over the remainder of 
the airfoil. Thus, the appropriate airfoil boundary condition in the direct 
region near the leading edge is the surface tangency requirement and in the 
inverse region, where the pressure is specified, it is essentially the specifi- 
cation of the derivative of the perturbation potential in the x-direction. In 
order to satisfy these at the airfoil boundary, which in general will not coin- 
cide with the Cartesian grid points, the derivatives in the boundary conditions 
are expanded as two term Taylor series about dummy points inside the airfoil. 

The derivatives in these series are then written in finite difference form 
using second order formulas for all first derivatives and at least first order 
ones for higher derivatives. In the direct region, central differences are 
used for x-derivatives and forward (on the upper surface) for the y-derivatives . 
However, to prevent numerical instability, the inverse region uses a second or- 
der backward difference formula for the first term of the Taylor series rep- 
resenting the x-derivative. 

In the inverse case the airfoil must also be computed by integrating the 
surface tangency condition for the ordinates, y, as a function of x, with the 
initial conditions given by the" slope and surface ordinate at the interface 
between the direct and inverse regions. In the present method, the latter are 
known because the nose region is solved directly, and the integration is accom- 
plished using the Runge-Kutta method of order four. 

It is believed, based on comparison with other results (ref. 3), that this 
method is an accurate and numerically consistent approach to the design and 


1388 



analysis of transonic airfoils in inviscid flow. For further details concerning 
the finite difference, boundary condition, formulations, etc., see references 
1 and 2. 


BOUNDARY LAYER ANALYSIS 

Experimental evidence (ref. 5-6) indicates that viscous boundary layer 
effects are very important in transonic flow. For example, the difference be- 
tween the actual airfoil surface and the effective surface, i.e. the displace- 
ment surface, can cause an airfoil inviscidly designed to have a lift coef- 
ficient of 0.6 to actually develop 25%-50% less lift. To prevent such discrep- 
ancies, the effect of the boundary layer displacement thickness, S*, should be 
included in both the analysis and design portions of any numerical method. 

In the present approach, the basic idea in the design case is to treat the 
airfoil determined by the inverse method as the displacement surface and to 
subtract from it the displacement thickness determined by a boundary layer com- 
putation. The result should be the actual airfoil ordinates. For the analysis 
case, the approach is to calculate a boundary layer displacement thickness and 
to use it to correct the location of the displacement surface (i.e. airfoil or- 
dinate plus 6*). The inviscid flowfield is then solved as before (ref. 1), 
where at present the displacement surface is updated every ten relaxation cycles. 

Obviously the boundary layer scheme must be efficient, reliable, and 
accurate. Thus, three integral methods were considered for inclusion in the 
present numerical method — Walz Method II (ref. 7), the Nash-Macdonald method 
with smoothing (ref. 8-9), and Green's lag-entrainment method (ref. 10). 

Figure 1 compares for a typical transonic case the upper surface displacement 
thickness predictions from these methods with those obtained by Bavitz 
(ref. 11), who used the Bradshaw scheme (ref. 12) modified with a trailing edge 
correction. (The Walz results are not plotted but are between the Bradshaw and 
Nash-Macdonald data.) Notice that the predictions are essentially identical 
over most of the airfoil, and all but Green's method predict separation near 
the trailing edge. Apparently Green's method needs some type of trailing edge 
correction. Since the Walz and Green methods numerically fail at separation 
due to their empirical equations and since the Nash-Macdonald approach is 3-6 
times faster, it was selected for incorporation into the present transonic air- 
foil design-analysis program. 


ANALYSIS RESULTS 

To start the direct problem the airfoil shape is inputed and a cubic 
spline fit and used to determine the ordinates and slopes in the computational 
plane. Next, to get some reasonable perturbation potentials, fifty relaxation 
cycles are performed on a very coarse grid (typically 13x7). Then the grid 
spacing is halved to a coarse grid (25 x 13) and fifty more cycles computed. 

At that point 6* is computed and the displacement surface ordinates updated 
using under relaxation, i.e. <S* n ew = 5* 0 ld + w (^* ~ fi*old) ■ The slopes are then 
determined from cubic splines through the new ordinates. The ordinates are 
updated every ten cycles thereafter. Typically, 400 cycles are performed on 
the coarse grid before halving to the medium grid (49 x 25) , where 200-250 
cycles are carried out. While this grid yields 66 points on the airfoil, the 


1389 



grid may be, if desired, halved again (to 97 x 49) to obtain 130 points on the 
airfoil. Compared to inviscid cases, convergence on the fine grid is slow, and 
400 relaxation cycles may be required. Fortunately, in many cases accurate 
results can be obtained on the medium grid. However, if double shocks exist, 
the fine grid may be needed to resolve them accurately. Convergence is deter- 
mined by monitoring the changes in perturbation potential and 6*. 

Typical total computation times on an Amdahl 470/V6 are one minute for 
medium grid results and less than 4 minutes for fine grid data (about 10 min- 
utes on a CDC 6600) . Convergence is usually faster on CDC type machines due to 
the increase in significant digits. 

As shown on figure 2, the present method can be used to demonstrate the 
effects of viscous interaction on a Korn 75-06-12 airfoil (ref. 9) near its 
design point. While the primary effect of the boundary : -yer is a 25% decrease 
in lift from the inviscid design value due to the pressure change on the upper 
surface, the lower surface pressure distribution is also affected.- In addition, 
notice the excellent agreement between the viscous theory predictions and NAE 
wind tunnel data (ref. 6). 

Another comparison with experimental data is shown on figure 3, again for 
the Korn airfoil but at an off-design condition. Although the minimum peak 
pressure and shock jump are slightly in error, which is not surprising since 
the method uses nonconservative finite differences, the overall agreement is 
excellent. Other results are shown in Table 1, which compares data obtained 
with the present theory with experimental values at about the same lift. Many 
of these were obtained using the medium grid results only, and it is believed 
that they indicate that the present viscous analysis method is adequate for 
engineering studies of lift, drag, and moment variation. This is particularly 
true when the large discrepancies between viscous and inviscid results and the 
small differences between viscous results and experiments are considered. 

DESIGN RESULTS 

The computational procedure ns£d in the inverse case is the same as de- 
scribed in references 2 and 3 except that 5* is subtracted from the displace- 
ment surface. The final design is usually obtained on the medium grid after 
250 relaxation cycles, although strong aft-cambered cases may require 400. 

At present there are some minor difficulties, due to the use of backward dif- 
ferences on the pressure boundary condition, in obtaining a desired pressure 
distribution near the trailing edge; and thus, the possibility of using the 
fine grid in design is under study. At present, a typical inverse run takes 
about 70 seconds (Amdahl 470/V6) . It should be noted that in obtaining the 
actual airfoil ordinates a transition point must be selected, and the final 
shape is somewhat sensitive to this choice. 

The importance of including the boundary layer in the design process is 
shown on figure 4. The shockless pressure distribution (solid line) used for 
this typical inverse design yielded for this Mach 0.72 case a 16% thick, highly 
aft-loaded airfoil. If the displacement surface ordinates were to be used to 
fabricate the airfoil instead of the correct values (i.e. if the boundary layer 
were ignored) , the actual pressure distribution would be as shown by the 
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symbols. This result, which was obtained using the present viscous analysis 
method, shows a lift and moment 20% less than the design values. 

On the other hand, when 6* is included in determining the ordinates and 
the resultant airfoil is analyzed with viscous effects included, the agreement 
is much better. The medium grid results for this case are shown on figure 5. 
Here, transition is assumed to occur just aft of the minimum peak pressure; and 
the 6* computed for the design pressure distribution differs by less than 

0.013% from the 6* determined in the analysis calculation. Nevertheless, as 
can be seen on figure 5, there is still a slight difference in the pressure 
distributions. However, considering the accuracy of boundary layer computa- 
tions, sensitivity to transition location, and the large difference between 
viscous and inviscid results, the agreement is quite good. In addition, it 
indicates acceptable numerical consistency between the present analysis and 
design techniques. 


CONCLUSION 

Based upon experimental comparisons, it is believed that the present 
viscous analysis method is suitable for obtaining engineering estimates of the 
characteristics of transonic airfoils. In addition, while the inverse design 
method has not been verified, it has been shown to be numerically consistent 
with the analysis results at the medium grid level. Efforts to extend the 
design procedure to a fine grid are in progress and will be reported later. 
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Table I.- COMPARISON WITH EXPERIMENT 




c 


C 



C 




D 


D 



D 



M 

CO 

Theory 

Expt 

Theory 

Expt 

Theory 


Expt 

1 

.622 

.48 

.49 

-.105 

-.104 

.0070 


.0070 

2 

.699 

.667 

.667 

-.107 

-.106 

.0108 


.0107 

3 

.75 

-.097 

-.093 

-.111 

-.114 

.0209 


.0146 

4 

.75 

.104 

.100 

-.118 

-.126 

.0086 


.0092 

5 

.75 

.454 

.458 

-.120 

-.121 

.0086 


.0090 

6 

.75 

.516 

.515 

-.122 

-.121 

.0089 


.0086 

7 

. 75 

.523 

.523 

-.119 

-.120 

.0092 


.0094 

8 

.75 

.539 

.535 

-.122 

-.120 

.0095 


.0088 

9 

.75 

.601 

.597 

-.119 

-. 122 

.0110 


.0081 

10 

.75 

.711 

.712 

-.130 

-.125 

.0151 


.0125 

11 

.15 

.475 

.47 

-.099 

-.105 

.0085 


.009 

Note: 

Cases 

1-11, 75-06-12 

at RN=2 1 

x 10 6 , Case 

11 GA(W)-2 

at 6 x 

10 6 
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Figure 1.- Comparison of boundary layer results for 
Korn airfoil 75-06-12. 
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